---
title: "The effect of radiation on the gut bacteriome of Aedes albopictus"
author: "Dongjing Zhang, Shi Chen, Adly Abd-Alla and Kostas Bourtzis"
date: "25/04/2021"
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## load data and needed packages

```{r}
setwd("C:/Users/abdallaa/OneDrive - IAEA/My/Kostas_paper 2020/Final_version_Feb_2021/Revission/raw_data")
alpha <- read.csv("Fig_1.csv")

str(alpha)
attach(alpha)
head(alpha)
alph=na.omit(alpha)
alph


library(ggplot2)
library(lattice)
library(gcookbook)
library(datasets)
library(MASS)
library(survival)
library(rmarkdown)
library(knitr)
library(coxme)
library(lme4)
library(nlme)
library(gapminder)
library(rcompanion)
library(FSA)
library(ggthemes) # Load
```

## Statiscic for the manuscript (in sequence)


```{r}
wilcox.test(Richness~Treatment1, data=alph)

wilcox.test(Richness~Sex1, data=alph)

kruskal.test(Richness~Age1, data = alph)
tmp <- dunnTest(Richness~Age1, data = alph, method ="none")
tmp


kruskal.test(Pielou.s.evenness~Sample, data = alph)
tmp <- dunnTest(Pielou.s.evenness~Sample, data = alph)
tmp

kruskal.test(Shannon~Sample, data = alph)
tmp <- dunnTest(Shannon~Sample, data = alph, method ="none")
tmp

kruskal.test(Simpson~Sample, data = alph)
tmp <- dunnTest(Simpson~Sample, data = alph)
tmp

kruskal.test(Shannon~Age1, data = alph)
tmp <- dunnTest(Shannon~Age1, data = alph, method ="none")
tmp

wilcox.test(Shannon~Treatment1, data=alph)

wilcox.test(Shannon~Sex1, data=alph)
```

## Statisctics for preparing table 1

```{r}
#Pielou.s.evenness
kruskal.test(Pielou.s.evenness~Sample, data = alph)
tmp <- dunnTest(Pielou.s.evenness~Sample, data = alph, method ="none")
tmp

#Richness
kruskal.test(Richness~Sample, data = alph)
tmp <- dunnTest(Richness~Sample, data = alph, method ="none")
tmp

#shannon
kruskal.test(Shannon~Sample, data = alph)
tmp <- dunnTest(Shannon~Sample, data = alph, method ="none")
tmp

#Simpson
kruskal.test(Simpson~Sample, data = alph)
tmp <- dunnTest(Simpson~Sample, data = alph, method ="none")
tmp
```

## To prepare Figure 1

```{r}
fig1a<-ggplot(alph, aes(x=Treatment1,y=Richness, fill=Treatment1)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2) + 
  ylim(20,60)
fig1a
tiff("fig1a.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(fig1a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = 'none') + xlab(expression(bold("Treatment")))+ ylab(expression(bold("Richness_index"))))
dev.off()

kruskal.test(Richness~Treatment, data = alph)

#Fig_1b

fig1b<-ggplot(alph, aes(x=Age1,y=Richness, fill=Age1)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2) + 
  ylim(20,60)
fig1b
tiff("fig1b.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(fig1b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = 'none') + xlab(expression(bold("Age")))+ ylab(expression(bold("Richness_index"))))
dev.off()

kruskal.test(Richness~Age, data = alph)
tmp <- dunnTest(Richness~Age, data = alph, method ="none")
tmp

#Fig_1c

fig1c<-ggplot(alph, aes(x=Sex1,y=Richness, fill=Sex1)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2) + 
  ylim(20,60)
fig1c
tiff("fig1c.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(fig1c+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = 'none') + xlab(expression(bold("Sex")))+ ylab(expression(bold("Richness_index"))))
dev.off()

kruskal.test(Richness~Sex, data = alph)

#Fig_1d

fig1d<-ggplot(alph, aes(x=Treatment1,y=Shannon, fill=Treatment1)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2) + 
  ylim(1,3.5)
fig1d
tiff("fig1d.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(fig1d+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = 'none') + xlab(expression(bold("Sex")))+ ylab(expression(bold("Shannon_index (Diversity)"))))
dev.off()

with(alph,boxplot(Shannon~Treatment))

#Fig_1e
fig1e<-ggplot(alph, aes(x=Age1,y=Shannon, fill=Age1)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2) + 
  ylim(1,3.5)
fig1e
tiff("fig1e.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(fig1e+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = 'none') + xlab(expression(bold("Age")))+ ylab(expression(bold("Shannon_index (Diversity)"))))
dev.off()

kruskal.test(Shannon~Age, data = alph)
tmp <- dunnTest(Shannon~Age, data = alph, method ="none")
tmp

#Fig_1f
fig1f<-ggplot(alph, aes(x=Sex1,y=Shannon, fill=Sex1)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2) + 
  ylim(1,3.5)
fig1f
tiff("fig1f.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(fig1f+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = 'none') + xlab(expression(bold("Sex")))+ ylab(expression(bold("Shannon_index (Diversity)"))))
dev.off()

kruskal.test(Shannon~Sex, data = alph)
```
## To prepare supplemetary table 3
## Compare the indixes between the two regions for all samples

```{r }
setwd("C:/Users/abdallaa/OneDrive - IAEA/My/Kostas_paper 2020/Final_version_Feb_2021/Revission/raw_data")
da1 <- read.csv("table_S3.csv")
#fig1$Time=as.factor(fig1$Time)
da1=na.omit(da1)
head(da1)

kruskal.test(Pielou.s.evenness~Regions, data=da1)
kruskal.test(Richness~Regions, data=da1)
kruskal.test(Shannon~Regions, data=da1)
kruskal.test(Simpson~Regions, data=da1)
```

## Compare the indicies between the two regions for each samples

```{r }
da2 <- subset(da1, Sex=="Male" & Age=="Pupa" & Treatment=="Irradiated")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da2)
kruskal.test(Richness~Regions, data=da2)
kruskal.test(Shannon~Regions, data=da2)
kruskal.test(Simpson~Regions, data=da2)

da3 <- subset(da1, Sex=="Male" & Age=="Pupa" & Treatment=="Control")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da3)
kruskal.test(Richness~Regions, data=da3)
kruskal.test(Shannon~Regions, data=da3)
kruskal.test(Simpson~Regions, data=da3)

#=============================================================
da4 <- subset(da1, Sex=="Female" & Age=="Pupa" & Treatment=="Irradiated")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da4)
kruskal.test(Richness~Regions, data=da4)
kruskal.test(Shannon~Regions, data=da4)
kruskal.test(Simpson~Regions, data=da4)

da5 <- subset(da1, Sex=="Female" & Age=="Pupa" & Treatment=="Control")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da5)
kruskal.test(Richness~Regions, data=da5)
kruskal.test(Shannon~Regions, data=da5)
kruskal.test(Simpson~Regions, data=da5)
#====================================================================

da6 <- subset(da1, Sex=="Male" & Age=="1D" & Treatment=="Irradiated")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da6)
kruskal.test(Richness~Regions, data=da6)
kruskal.test(Shannon~Regions, data=da6)
kruskal.test(Simpson~Regions, data=da6)

da7 <- subset(da1, Sex=="Male" & Age=="1D" & Treatment=="Control")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da7)
kruskal.test(Richness~Regions, data=da7)
kruskal.test(Shannon~Regions, data=da7)
kruskal.test(Simpson~Regions, data=da7)

#=============================================================
da8 <- subset(da1, Sex=="Female" & Age=="1D" & Treatment=="Irradiated")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da8)
kruskal.test(Richness~Regions, data=da8)
kruskal.test(Shannon~Regions, data=da8)
kruskal.test(Simpson~Regions, data=da8)

da9 <- subset(da1, Sex=="Female" & Age=="1D" & Treatment=="Control")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da9)
kruskal.test(Richness~Regions, data=da9)
kruskal.test(Shannon~Regions, data=da9)
kruskal.test(Simpson~Regions, data=da9)
#====================================================================

da10 <- subset(da1, Sex=="Male" & Age=="4D" & Treatment=="Irradiated")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da10)
kruskal.test(Richness~Regions, data=da10)
kruskal.test(Shannon~Regions, data=da10)
kruskal.test(Simpson~Regions, data=da10)

da11 <- subset(da1, Sex=="Male" & Age=="4D" & Treatment=="Control")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da11)
kruskal.test(Richness~Regions, data=da11)
kruskal.test(Shannon~Regions, data=da11)
kruskal.test(Simpson~Regions, data=da11)

#=============================================================
da12 <- subset(da1, Sex=="Female" & Age=="4D" & Treatment=="Irradiated")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da12)
kruskal.test(Richness~Regions, data=da12)
kruskal.test(Shannon~Regions, data=da12)
kruskal.test(Simpson~Regions, data=da12)

da13 <- subset(da1, Sex=="Female" & Age=="4D" & Treatment=="Control")
head(da2)

# kruscal.tset
kruskal.test(Pielou.s.evenness~Regions, data=da13)
kruskal.test(Richness~Regions, data=da13)
kruskal.test(Shannon~Regions, data=da13)
kruskal.test(Simpson~Regions, data=da13)
```
